Non-conformally flat initial data for binary compact objects 
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A new method is described for constructing initial data for a binary neutron-star (BNS) system 
in quasi-equilibrium circular orbit. Two formulations for non-conformally flat data, waveless (WL) 
and near-zone helically symmetric (NHS), are introduced; in each formulation, the Einstein-Euler 
system, written in 3+1 form on an asymptotically flat spacelike hypersurface, is exactly solved for 
all metric components, including the spatially non-conformally flat potentials, and for irrotational 
flow. A numerical method applicable to both formulations is explained with an emphasis on the 
imposition of a spatial gauge condition. Results are shown for solution sequences of irrotational BNS 
with matter approximated by parametrized equations of state that use a few segments of polytropic 
equations of state. The binding energy and total angular momentum of solution sequences computed 
within the conformally flat — Isenberg- Wilson-Mathews (IWM) — formulation are closer to those of 
the third post-Newtonian (3PN) two point particles up to the closest orbits, for the more compact 
stars, whereas sequences resulting from the WL/NHS formulations deviate from the 3PN curve 
even more for the sequences with larger compactness. We think it likely that this correction reflects 
an overestimation in the IWM formulation as well as in the 3PN formula, by ~ 1 cycle in the 
gravitational wave phase during the last several orbits. The work suggests that imposing spatial 
conformal flatness results in an underestimate of the quadrupole deformation of the components of 
binary neutron-star systems in the last few orbits prior to merger. 
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I. INTRODUCTION 

Inspiral to merger of binary neutron stars (BNS) 
is one of the most promising sources of ground-based 
gravitational-wave detectors. A fully general relativis- 
tic numerical simulation is the unique approach to pre- 
dict the gravitational waveform from the late inspiral to 
merger phase. Such a simulation begins with preparing 
quasi-equilibrium initial data with a close orbital separa- 
tion ~ 45-50 km. 

Quasi-equilibrium initial data for binary neutron stars 
introduce two kinds of inaccuracies into inspiral simu- 
lations. One is due to ignoring the radial component 
of the velocity of orbiting stars, the other to artificial 
restrictions on the geometry of the initial hypersurface 
A common choice for the geometry of the initial 
hypersurface is a conformally flat three-geometry [EH, 
and a similarly restrictive alternative is presented in [J]. 1 
The former error is reduced by adding radial velocity to 
minimize the oscillation around the inspiral orbit, where 
the radial velocity may be determined empirically or cal- 
culated from the post-Newtonian formula of inspiraling 
point masses. Both errors become negligible if the initial 
separation of the binary is large enough, possibly five or- 
bits or more before the merger; but increasing separation 



increases the cost of computing time, and maintaining 
accuracy in numerical simulations may still be an issue. 2 

In a previous letter Q, we have reported that the in- 
accuracy of the binary orbit arising from spatial confor- 
mal flatness can be largely removed if one solves the full 
Einstein equation for all metric components, including 
the non-conformally flat part of the spatial metric, on 
a Cauchy surface Et, using the formulation presented in 
H> B 03 • In this formulation, Einstein's equation is writ- 
ten in a 3+1 form and the time derivative of the confor- 
mal three metric, dt^ab, which carries the dynamics of 
the spacetime in our choice of the gauge, is set to zero. 
7 Q ;, is conformally related to the spatial metric ^ a b in each 
slice S t by 7 a b = ip^Jab, with tp a conformal factor. As 
a result, the field equations for the metric components 
become elliptic equations on an initial slice S t , and they 
yield an asymptotically flat metric. We call this approach 
the waveless formulation (WL). 

We have also experimented with another formulation 
for quasi-equilibrium initial data in which all components 
of the metric are computed; preliminary results were pre- 
sented in [l3T |. In this approach, helical symmetry is im- 
posed in the near zone from the center of mass to the 
radius ~ A = n/Q, and either the WL formulation is 
applied outside, or the computational domain is trun- 



1 For the computation of black hole-neutron star binary in quasi- 
equilibrium, see e.g. Q 



2 For long-term simulation of binary black hole inspirals and 
matching to the post-Newtonian results, see e.g. [6( 
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cated at this radius. Here, fl is the orbital angular veloc- 
ity and A is the wavelength of the dominant, primarily 
£ = m = 2 quadrupole, mode of the gravitational waves 
expected to be radiated from the system. In this paper, 
we discuss the near-zone helically symmetric (NHS) for- 
mulation together with the WL formulation, for which 
numerical methods are common. 

A significant difference in the binding energy and total 
angular momentum between the solutions from WL /NHS 
formulations and those from a conformally flat formula- 
tion, the Isenberg- Wilson-Mathews (IWM) formulation 
[TTI . Il2j , is found and is discussed in the later section. In 
the IWM formulation, the solution sequences - the plots 
of these quantities as functions of orbital angular velocity 
- become closer to those of third post-Newtonian (3PN) 
point particles up to the closest orbits when the com- 
pactness of each star is increased. In contrast, sequences 
obtained from the WL/NHS formulations deviate more 
from the 3PN curve for larger compactness. 

We expect waveless and helically symmetric solutions 
to accurately approximate the outgoing metric in the 
near zone, where the gravitational wave amplitude is 
small compared to the Coulomb part of each metric po- 
tential. Results of Ref. [13[ and of the present paper 
support this expectation by showing that corresponding 
WL and NHS solutions nearly coincide. 

Several groups have developed simulation codes for 
BNS inspirals and merger; stable long-term simulations 
[LH [l5j . magnetized BNS simulations 0], and black 
hole-neutron star binary merger simulations (l7j are now 
feasible. As mentioned above, however, accurate model- 
ing of the last several orbits of inspiraling binary compact 
object using quasi-equilibrium sequences will be still use- 
ful, because the lower computational cost allows one to 
study gravitational wave sources by exploring a wider pa- 
rameter space, varying the mass ratio and the dense mat- 
ter EOS. One of the other applications will be the com- 
parison with the results of simulations, which becomes a 
reliable calibration for both of the numerical solutions. 

This paper is organized as follows: In Sec|TT] we de- 
scribe the WL/NHS formulations. These are essen- 
tially identical to those introduced in our previous pa- 
pers 0, [l(| EH , except for a few modifications suitable 
for coding. All equations used in actual coding are writ- 
ten in Appendix [Al in detail. As a model for the EOS of 
high density matter, the parametrized EOS developed in 
[HI [H[ is used in the computations and is briefly intro- 
duced in this section. In Sec lIIII the numerical method is 
discussed, with emphasis on the major differences from 
the previous conformally flat code. In Sec lIVI we report 
results from the WL/NHS computation of binary systems 
and of constant rest mass quasi-equilibrium sequences. 

In this paper, spacetime indices are Greek, spatial in- 
dices Latin, and the metric signature is — I — I — h For writ- 
ing the basic equations, geometric units with G = c = 1 
are used, while for tabulating the numerical solutions, 
cgs units or other appropriate units are used. 



II. FORMULATION 

A. 3+1 decomposition and gauge conditions 

The spacetime M. = K x £ is foliated by a family of 
spacelike hypersurfaces (T, t )tem parametrized by t. The 
future-pointing unit normal n a to the hypersurface St is 
related to the generator t a of time translations, for which 
t a V a t = 1, by 

t a ^an a +(3 a . (1) 

Here a is the lapse function and j3 a the shift vector, with 
(3 a n a = 0. n a is related to the gradient of t by n a — 
—aV a t. It is also related to the helical vector k a , the 
generator of time translation in a rotating frame, by 

k a =an a +uj a , (2) 

where a spatial vector oj a := (3 a + Q,(j) a is the rotating 
shift in the rotating frame, and £1 is a constant angular 
velocity of the rotating frame. The helical vector k a is 
not everywhere timelike, but it is transverse to the sur- 
face E t , and normalized as k a Vt = 1. 

The spatial metric iab(t) induced on E t by the space- 
time metric g a p is equal to the projection tensor orthog- 
onal to n Q , 7 Q /3 = g a f3 + n a np, restricted to E t . We 
introduce a conformal factor tp, a conformally rescaled 
spatial metric j a b, and a flat spatial metric f a b, with 
lab = ip^lab, and with the conformal factor specified by 
the condition 7 = /, where 7 and / are the determinants 
of jab and f a b. In a chart (t,x a ), the metric g a p has the 
form 

ds 2 = -a 2 dt 2 + ip% b (dx a + (3 a dt)(dx b + (3 b dt). (3) 

Let us denote by h a b and h ab the differences between the 
conformal metric and the flat one: 

lab = fab + h ab , I 11 " - + h ab . (4) 

The extrinsic curvature of each slice S t is defined by 

K a b = -^nlab, (5) 

where the action of £ n on j a b in the above definition, 
and on other spatial tensors hereafter, is given by 

£n7afc := —dtlab ~ ~^plab\ (6) 

a a 

here d t i a b is the pullback of £tlap to S t , with £ t the Lie 
derivative along the vector t a defined on A4, and £p is 
the Lie derivative along the spatial vector (3 a on S t . 
Einstein's equation is written in the 3+1 form 



[G aP - 8TrT aP )n a n' 3 = 0, (7) 

(G Q/3 - 8TrT af3 )i a a nP = 0, (8) 

(G a0 - 8irT af3 ) (V + \n a n^ = 0, (9) 

(Gap - 8irT af3 ) ( 7^ " llabl al3 ) = 0. (10) 
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These equations are the Hamiltonian and momentum 
constraints, the trace of the spatial projection combined 
with the Hamiltonian constraint, and the tracefree part 
of the spatial projection, respectively. They are solved 
for tp, [3 a , the combination aifi, and h a b- For perfect-fluid 
spacetimes, the stress-energy tensor T a ^ is written 



T afi = (e + p)u a u f:i +pg al3 , 



(11) 



where e is the energy density, p the pressure, and u a the 
4-velocity of the fluid. 

The above set of equations are solved imposing as co- 
ordinate conditions the maximal slicing condition, 

K = 0, (12) 

and the generalized Dirac gauge condition [S [IS EH , 

D b r" = D b h ab = 0, (13) 

o 

where D a is the covariant derivative associated with the 
fiat metric f a b- Concrete forms of Eqs.(j?j)- (fTU)) are pre- 
sented in Appendix |A"1 



B. Waveless and near-zone helically symmetric 
formulations 

As a model for binary compact objects in general rel- 
ativity, helically symmetric spacetimes have been intro- 
duced [HI and studied by several authors [IS HH US US 
US US US H3- Helically symmetric binary solutions for 
point-particles in a post-Minkowski framework 28] anal- 
ogous to the electromagnetic two-body solution [29| , and 
for several toy models have been calculated [l3l |. 

Helically symmetric spacetimes do not admit fiat 
asymptotics. However, it is expected that, up to a certain 
truncation radius where the energy of radiation does not 
dominate the gravitational mass of the system, solutions 
have an approximate asymptotic region in which gravi- 
tational waves are propagating in a curved background. 
Such a solution, however, has not yet been calculated 
successfully in the regime of strong gravity. 

Helical symmetry, 



(14) 



implies for the 3-metric and extrinsic curvature on a ini- 
tial hypersurface E t , 

£fc7ab = 0, £kK a b = 0. (15) 
Using the relation k a = an a + uj a , we have 



a 



£ K 



-"n-^ab -^cu^ ab- 
Ct 



(16) 
(17) 



Because k a is timelike in the fluid, helical symmetry for 
the fluid variables, 



has the meaning of stationarity for a rotating observer. 

Our formulation for the non-conformally flat data of 
binary compact objects in a quasi-equilibrium quasi- 
circular orbit is based on the helically symmetric formu- 
lation. We further impose either a waveless condition or 
near-zone helical symmetry in the gauge (TT2")) and (Tj"3"|) . 

a. Waveless formulation As discussed in [lOj], the 
condition, dt x f ab — 0(r~ 3 ), is sufficient to enforce 
Coulomb- type fall off in the asymptotics. For our wave- 
less formulation in this paper, we impose the stronger 
condition 

dtlab = 0, (19) 
which amounts to writing the extrinsic curvature as 



K„ 



1 

2^ 
1 

2^" 



£/37ab 
£/37ab 



1 

2~T\ 
1 

2^" 



lab 



jab Q£<t> Ini/' 4 



(20) 



where helical symmetry is used to get the second equality. 
Only the first term on the r.h.s. remains in the maximal 
slicing condition. Because the trace of Eq. (f2"0"|) has the 
same form for K = as the trace of the original equation 
(J5|), the waveless condition (fT9|) does not affect the max- 
imal slicing condition. Note that the second term of the 
r.h.s. of Eq. ([2U)l does not appear in the tracefree part of 
K a b; in other words, the time derivative of the conformal 
factor ip does not appear in the initial value formulation 
in this slicing. The other time derivatives are given by 
the helical symmetry conditions, Eqs. tfTT)) and (fTS]) . 

b. Near-zone helically symmetric formulation Near- 
zone helical symmetry means that we impose helically 
symmetric conditions (|16 j) -(|18 |) in the region from the 
center of the source to about one wavelength of the I = 
m = 2 mode of the gravity wave, r < A :— tt/SI; we then 
either truncate the domain of numerical computation at 
this radius or use the waveless formulation outside. The 
latter implies for K a b the condition 



K„ 



la 



for r < aX, 



77-£/37q& + ^— Ta&^JE^lnV" 4 for r > aA, 
2a 2a 



(21) 



£ k u a = 0, £ fc e = 0, £ fc ;p = 0, 



(18) 



where the constant a, the coordinate radius of the heli- 
cally symmetric zone in units of A = n/il, is restricted 
to a < 1.5. Without this restriction, iterations fail to 
converge to a binary solution. In the near-zone-helical 
+ outside- waveless formulation, all metric components, 
including those of the spatial metric, have Coulomb-type 
fall off. We have compared the NHS solution to the WL 
solution in our previous paper [l3l ] and confirmed that 
the difference in the non-conformal flat part of the spa- 
tial metric is about 1% for the BNS of M\/R ~ 0.17, 
where Mi/R is the compactness, the ratio of the gravi- 
tational mass to the circumferential radius of a spherical 
star having the same rest mass as each component star 
of the binary. 
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C. Formulation for the irrotational flow 

The late stage of BNS inspiral is modeled by a con- 
stant rest mass sequence of quasi-equilibrium solutions 
with negligible spins and magnetic fields, a description 
appropriate to a binary of old pulsars with spin periods 
longer than 100 ms. Since the viscosity of the high den- 
sity matter is expected to be negligible, a neutron star in 
a binary system is not spun up by the tidal torque dur- 
ing the inspirals. Hence, the flow field remains approxi- 
mately irrotational, and each neutron star is modeled by 
an irrotational perfect fluid [30j. 

The equation of motion, V pT a P = 0, for a perfect fluid 
has the form, 



P 



hu a Vf3(pu p ) - pTV a s = 0, (22) 



where s is the entropy per baryon mass, h is the relativis- 
tic enthalpy per baryon mass h := (e + p)/p, and local 
thermodynamic equilibrium dh = Tds + dp/ p is assumed. 
We assume constant entropy per baryon (s = const) ev- 
erywhere inside the neutron star, together with a one- 
parameter EOS, 



The form 



p=p(p). 



u Vp(hu a ) + W a h = 



(23) 



(24) 



of the relativistic Euler equation then follows from local 
conservation of baryon mass, 



V a (pu a ) = 0. 



(25) 



Written in terms of the Lie derivative along u a , these last 
equations have the form 



-±=£ u {p^g) = 0, 
£ u {hu a ) + V a h = 0. 



(26) 
(27) 



A state is stationary state in the rotating frame if it is 
helically symmetric, if each physical field is Lie derived 
by the helical vector field k a , as in Eq. (jTHJ), or 



£fe(pit*y'— g) = 0, and £k(hu a ) = 0, 



(28) 



where w* is the scalar u a V a t. 

The relativistic Euler equation (f24|) can be rewritten 

as 



0, 



where 



up a := Vp(hu a ) - V a {hufs) 



(29) 



(30) 



is the relativistic vorticity tensor. This implies that, for 
irrotational flow, hu a has a potential $, 



htlr 



(31) 



and hence the relativistic Euler equation has a first in- 
tegral. With a spatial velocity v a in the rotating frame 
defined by 



u a = u\k a +v a ), 
where v a n a = 0, Eq. (f2"T)) becomes, 



(32) 



£ u (hu a )+X7 a h = u 



£ k + v (hu a ) +V a — 



w'V^ £„$ 



therefore the first integral is 

+ - 



0; (33) 



(34) 



where £ is a constant 3 . Note that Eqs. (f2"8"| and (|3"Tj) 
imply a flow with £fc<!> = constant. Such a flow is both 
irrotational and helically symmetric with the shape of 
the star fixed in the rotating frame. Solutions describing 
irrotational binaries in Newtonian and post-Newtonian 
gravity are found in (3lj . and details of the formulation 
for helically symmetric irrotational flow are given in [20L 

There are three fluid variables and two parameters 
to be determined in the above formulation. The fluid 
variables are a thermodynamic variable, the velocity po- 
tential $, and the time component of the 4- velocity u*; 
and these are calculated from the first integral (|34[) . 
the rest-mass conservation equation (|26p . and the nor- 
malization of the 4- velocity, u a u a — — 1. A concrete 
form of these equations are presented in Appendix IA 31 
For the independent thermodynamic variable, we choose 
q := p/p, and other thermodynamic variables are deter- 
mined from the thermodynamic relations and the one- 
parameter EOS, which are briefly explained in the next 
section. The number of fluid variables and parameters 
are augmented in the numerical computation, which is 
mentioned in Appendix IB 21 (or see [371]). 



D. Parametrized equations of state 

Recently, a parametrization for the EOS of nuclear 
matter has been studied, and it is shown that a 
parametrized EOS with three polytropic intervals ap- 
proximates with fair accuracy a variety of current can- 
didate EOS, over a range of densities that extends from 



3 Cartan identity k@u>p a = Lj s (hu a ) — 'V a (huak^) implies, for the 
helically symmetric irrotational flow satisfying £k(hu a ) = and 
w /3a = 0, a relation, hu a k a = constant, equivalent to Eq. ( 1 3 4 1 ) . 
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the inner crust to the maximum neutron-star density [1 81 ] . 
Two of these intervals and three parameters cover den- 
sities below the central density of a 1.4M© neutron star, 
and waveforms from binary inspiral can be used to con- 
strain this three-dimensional subspace of the parameter 
space [l5j |. This parametrized EOS is used in our models 
for BNS data. 



1. Construction of piecewise polytropic EOS 

In presenting these piecewise polytropes, it is helpful 
to introduce a relativistic Emden function q by 



q :=p/p 



(35) 



and to write the remaining thermodynamic variables in 
terms of q. For an isentropic flow with s = 0, the local 
first law of thermodynamic equilibrium, dh = Tds + j^dp, 
takes the form 



dh 



-dp, 



where h is the enthalpy per baryon mass. 
A piecewise polytropic EOS is given by 



p = Kip 



r, 



(36) 



(37) 



in the intervals p 6 [pi-i, Pi), i = 1, - ■ ■ , N, with po = 
and pn — oo. In this section the subscript i denotes the 
ith interval, associated with a set of constants {Ti,Ki} 
with i = 1, • • • ,N, and labels the value of quantities at 
the higher density side of each interval, \pi—i,Pi). Be- 
cause we consider only continuous EOS, Pi,hi,€i and g» 
are the values of each of these quantities at density pi. 

The constant indices I\ arc N model parameters, and 
values of one thermodynamic variable at interfaces com- 
prise a set of N — 1 model parameters. A requirement 
that the pressure at the interface is continuous, 



r, 



K i+ ip 



(38) 



uniquely specifies values of K{ up to one free parameter, 
one of Ki of a specific ith interval, which is usually spec- 
ified by prescribing the value of pressure pi at the cor- 
responding interface density pi. Therefore we have 2N 
parameters for a parametrized EOS with N intervals. 

To compute other thermodynamic quantities from q, 
we use the following relations, valid in the ith interval, 



h - hi-i 



P = K i 1 

-1 

p = K i q 

r 



r«-i 

ph - p, 



{q - qi-i), 



(39) 
(40) 
(41) 
(42) 



where Eq. (|4"Tj) is obtained by integrating the relation 

(43) 

in the ith interval q £ fe-i, <ft]. Here, 



dh = —dp = — dq 

p Ti - 1 



r 



with ho — 1 and qo = 0. 



(44) 



2. Choice for the parameters 

In the latter sections, we present the results of quasi- 
equilibrium BNS solutions calculated using two types of 
parametrized EOS. The first EOS contains one free pa- 
rameter, which is used to estimate the accuracy of the 
measurement of the EOS parameter, and the neutron- 
star radius, by gravitational-wave observations of the 
inspirals of BNS [l5j]. The second EOS is a four- 
parameter fit to the candidates of neutron-star EOS. 
Those candidate EOS are tabulated nuclear EOS, and the 
parametrized EOS with four parameters approximates 
each candidate within the rms residual typically in the 
order of ~ 0.1%, and ~ 4.3% for the worst case |18l |. 

The parametrized EOS with one parameter uses two 
polytropic intervals. The lower density interval approx- 
imates the known subnuclear density EOS, the fixed 
crust EOS, around 0.1p nuc ~ p nuc by setting (T Q ,K Q ) = 
(1.35692,3.59389 x 10 13 ). Here, p DUC is the nuclear satu- 
ration density, and the constant Kq is in cgs units which 
give the pressure p in dyn/cm 2 . For the second poly- 
tropic interval at the higher density side, the adiabatic 
index is set Ti = 3. Then, the pressure p\ at the density 
pi = 10 14 ' 7 g/cm 3 is chosen as a parameter, and the di- 
viding density at the fixed crust and the next polytropic 
piece po is determined as the intersection of the two in- 
tervals. Further details are found in [To| . 

The four-parameter fit uses the same crust EOS as 
above, and three other polytropic intervals. The adia- 
batic indices of higher polytropic intervals {Fi,r2,r3}, 
and the pressure p\ at the interface between i — 1 and 2 
are chosen as fitting parameters, while the dividing den- 
sity po is evaluated in the same way as above, and other 
dividing densities are fixed as p\ — 10 14 ' 7 g/cm 3 and 
p 2 = 10 15 g/cm 3 . The EOS parameters and correspond- 
ing data for the spherical solutions are summarized in the 
later section IIV1 Further details for the four-parameter 
fit are found in [lSl ] . 



III. COMPUTATION 

A system of elliptic equations and algebraic relations 
are solved applying a self-consistent field iteration scheme 
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|34j . Recently, the convergence of such scheme for Newto- 
nian barotropic stars has been mathematically analyzed 
in [H. 

The WL/NHS code for the irrotational BNS presented 
in this paper is developed on top of the former BNS code 
in which the IWM formulation is used . Another ver- 
sion of the WL/NHS code, based on the triaxially de- 
formed rotating neutron-star code described in [37j |. has 
been developed, and its results are presented elsewhere. 
The numerical method used in these codes is briefly re- 
peated in Appendix [Bl 



A. Imposition of Dirac gauge 

The primary difference between the WL/NHS code 
and an IWM code is the computation of the non- 
conformally flat part of the spatial metric j ab = f a b + h ab . 
The conformal spatial metric j ab has to satisfy two con- 

ditions, 7 = / and — 0, which turn out not to be 

automatically satisfied when the spatial tracefree part 
of Einstein's equation (TIT)]) or its concrete form in the 
code, either Eq. (|A57|) or (|A60[) , is solved for h ab . To 
impose these conditions on "/ ab accurately, we first make 

a gauge transformation of h a b to satisfy D b h = 0, and 
we then correct the conformal factor to enforce the rela- 
tion 7 = tp 12 f a t each iteration cycle. Note that these 
two conditions are not explicitly imposed in Eq. (|A57[) or 
(IA60P , and they are violated mainly due to the numerical 
error of finite differencing. 

The gauge vector is calculated numerically by the fol- 
lowing procedure: a perturbation of the spatial metric 

hab 

o o 

hab -> hab ~ D a ^ b - D b £ a , (45) 

implies, to the same order, that the conformally rescaled 
metric with 7 = / satisfies 

hab hab - DaZb - D b Za + ^/okA-f . (46) 

We adjust h ab to this order to satisfy the Dirac gauge 
condition; namely, writing 

h' ab = h ab - D a & - D b £ a + ^f ab D c $ c , (47) 

we let h' ab satisfy the Dirac gauge condition to linear 

o 

order in h ab , D b h' ab — 0, which leads to 

A^a + ^D a D% = D b h ab . (48) 

This equation is solved by introducing the decomposi- 
tion 

ia = G a - jD a B, (49) 



which results in a set of elliptic equations, 

AG a =D b h ab , and AB = D a G a . (50) 

These equations (|50|) are solved using the same Poisson 
solver described in Appendix [B] and a solution is substi- 
tuted in Eq. (gSJ) and then in Eq. ijlTjl. In the r.h.s. of 
Eq. (|47p , h ab is calculated from the tracefree part of Ein- 
stein's equation, either Eq. (|A57|) or (|A60| . and it is 
replaced by h' ab , which satisfies the Dirac gauge condi- 
tion more accurately. We have also experimented with 
a transformation of the contravariant components of h ab 

o 

analogous to Eq. (|47|) . and let D b h ab — be satisfied; 
however, the results did not change. 

After the above gauge transformation, the condition 
7 = / is imposed by adjusting the conformal factor -0 to 

$ = *l> (j) 12 , (51) 

where 7' is the determinant of j' ab = f ab + h' ab . Note 
that, to impose 7 = /, we do not change the value of 
h' ab . These two corrections to h ab and ip are made once 
per iteration 

The other parts of the method of computation, includ- 
ing the iteration scheme, are common to our previous 
codes [H, H3, Hl| , which are briefly reviewed in Appendix 

m 

B. Coordinate and grid parameters 

The WL/NHS code uses two coordinate patches: a 
spherical patch, called the central coordinate system, 
on which the metric components are calculated, and a 
surface-fitted spherical coordinate patch on which the 
fluid variables are computed. The origin of the central co- 
ordinates (r, 8, (j)) is the mass center of the binary system, 
and that of the surface-fitted coordinates {rf,9f,<j>f) is 
the geometric center of the component star, where f/ is 
related to the radial coordinate 77 by f/ = rf/R(6f,<j>f) 
and R(0f,cj)f) is the surface of the star. We match the 
radial coordinate lines at (8, <fi) — (ir/2,0) of the central 
coordinates, and that of the surface-fitted coordinates 
(6 f ,cj) f ) = (tt/2,0), and set the 6 = and 6 f = lines 
to be parallel. Only the octant of the whole space for 
the central coordinate is solved, while a quarter for the 
surface-fitted coordinate. The spherical coordinates cor- 
respond to the Cartesian coordinates in the usual way; 
the (0,(j>) = (tt/2,0) line to the (positive) x-axis, the 
(7r/2,7r/2) line to the y-axis, and the 8 = line to the 
z-axis. 

The accuracy of the numerical solutions depends on the 
resolution of the finite differencing determined by the grid 
spacings (Ar, A8, A(f>), and the order of the truncation of 
multipole expansion ^ max . The latter is constrained by 
the resolution since the multipoles involved in the Green's 
function, which oscillates rapidly for the larger £, should 
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JVr 


Number of intervals An 


in 


r e [0,r b ] (CC). 


n r 


Number of intervals Ari 


in r e [0,r c ] (CC). 


N e 


Number of intervals A6i 


in 


9 G [0,71-/2] (CC). 


N^ 


Number of intervals A(j>i 


in 


4>€ [0.7T/2] (CC). 


N f . 


Number of intervals Afj 


in 


f f e [o,i] (SFC). 


Nj 


Number of intervals A6i 


in 


e f e [0,7r/2] (SFC). 


Nl 


Number of intervals A(j>i 


in 


<t>, e [o,7r] (SFC). 



TABLE I: Summary of grid parameters. (CC) stands for the 
central coordinates, and (SFC) for the surface-fitted coordi- 
nates. 



r b r c N r n r Ng 




Nf 


AT/ 


N f l f 

<p max 


W 4 R 5R 250 160 64 


64 40 


32 


32 


24 8 



TABLE II: Coordinate parameters, and the number of grid 
points used in this paper. Ro is the geometrical radius of the 
neutron star along the (6f,(f>f) = (tt/2, 0) line. Z max and Z^ax 
are the highest multipoles included in the Legendre expansion 
in the central and surface-fitted coordinates, respectively. 



be resolved on the grids. The radial grid spacing Ar of 
the central coordinates is equidistant for r € [0, r c ], and 
increases in geometric progression for r € [r c ,rb]. The 
grid spacings of the other coordinates are equidistant. 
For further details, see [H, H3, HI] . For the grid param- 
eters, we choose values listed in Table [TTJ Typically, 1 
cycle of iteration takes about 70 s for this grid setup us- 
ing a single core of Intel Xeon CPU X5450 with 3.00GHz 
clock. 



IV. QUASI-EQUILIBRIUM SOLUTIONS 



10 r i 10 rr 




-10 -5 5 10 -10 -5 5 10 



x/R x/flo 

FIG. 1: Contours of (h xx — h yy )/2 (left panels) and h zz (right 
panels) in the xy-plane. Top panels are those of a WL solu- 
tion with the EOS parameter HB and the orbital radius d/Ro 
= 1.5, where _Ro is the coordinate radius (a half of the diame- 
ter) of the neutron star along the :r-axis. Contours are drawn 
every 0.001 step, where the solid (dashed) contours in the top 
panels corresponds to positive (negative) values of hij . Thick 
dotted circles are the surface of neutron stars. Bottom panels 
are the contours of the 2PN asymptotic formula (|52|) calcu- 
lated for the two point masses assuming the same coordinate 
separation d/Ro — 1.5, and mass M — 1.35Mo x 2. Contours 
are also drawn every 0.001 step. In the left top and bottom 
panels, the outermost contour (interrupted by the boundary 
of the figure) corresponds to —0.002. In the bottom two pan- 
els, contours out of range (< —0.01 for the left panel, and 
^ ±0.005 for the right) are truncated. 



A. Behavior of hij for selected solutions 

Quasi-equilibrium solutions of irrotational BNS are 
calculated for the various sets of EOS parameters summa- 
rized in Table HTT1 As an example of the WL solutions, 
we present in Fig[T] contours of selected components of 
for the parametrized EOS HB, with orbital radius 
d/Ro = 1.5, where Ro is the coordinate radius (half the 
diameter) of the neutron star along the x-axis. For a 
qualitative comparison, contours are also shown for the 
leading order terms 0(r _1 ) of the asymptotic solution 
of h^ in a second order post-Newtonian (2PN) approx- 
imation with maximal slicing and a transverse-traceless 
gauge for hij, as derived in [39| (see, Eq. (5.30)), namely 

hij = ^ ji^ij + 2 nk ( nl !kj + n J I ki ) - ^n 4 n J / fefe 

3 1 5 1 

+ -n t n : >n k n l I k i + -Sijl kk - -5 lj n k n l I k i \ 

+0(r- 2 ) (52) 



where 

Iij — J px l x : 'd' i x, and ri 1 — — . (53) 

In the quadrupole integrals Iij , we substituted two 
1.35Mq point masses, separated by the same coordinate 
length as the above WL solution. The region shown in 
these figures does not extend far enough to have asymp- 
totic behavior, though the contours qualitatively agree. 

In Fig. [51 selected components of hij are plotted along 
the x-axis for the cases with parametrized EOS 2H, HB, 
and 2B, from top to the bottom panels. In each case the 
the orbital radius is again d/Ro = 1.5. In our models, the 
gravitational mass of the corresponding spherical star is 
Mi = 1.35M and M\/R of each EOS increases in the 
order of 2H, HB, 2B (see Table [Till), which is reflected 
by the increasing amplitude of hij . Here and after, the 
compactness of each component star in the binary system 
means the value of Mi /R for a single spherical star with 
the same rest mass. 

In the right panels, corresponding to the left panels, 
log-log plots of the h yy and h zz components are shown 
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Model 


log(po) 


log(pi) 


Ti 


r 2 


r 3 


Mi [M ] 


Mo [M @ ] 


R [km] 


Mi/R 


log(pc) 


2H 


13.847 


34.90 


3 


3 


3 


1.35 


1.4549 


15.224 


0.13097 


14.573 


HB 


14.151 


34.40 


3 


3 


3 


1.35 


1.4927 


11.606 


0.17181 


14.918 


2B 


14.334 


34.10 


3 


3 


3 


1.35 


1.5251 


9.7268 


0.20500 


15.141 


SLy 


14.165 


34.384 


3.005 


2.988 


2.851 


1.35 


1.4947 


11.469 


0.17385 


14.934 


APR1 


14.294 


33.943 


2.442 


3.256 


2.908 


1.35 


1.5388 


9.1385 


0.21819 


15.221 


FPS 


14.220 


34.283 


2.985 


2.863 


2.600 


1.35 


1.5055 


10.702 


0.18631 


15.038 


BGN1H1 


14.110 


34.623 


3.258 


1.472 


2.464 


1.35 


1.4789 


12.626 


0.15792 


14.912 


ALF3 


14.188 


34.283 


2.883 


2.653 


1.952 


1.35 


1.5069 


10.350 


0.19264 


15.150 



TABLE III: Parameters of each EOS and properties of the spherical neutron-star model based on that EOS and having 
gravitational mass Mi = 1.35M0. The pressure pi [dyn/cm 2 ] is the value at the dividing density pi = 10 14 ' 7 g/cm 3 , and values 
of log(pi) and {Pi, P2, T^} are taken from Table I of [15| and Table III of [18^. The parameters to fit the crust EOS are chosen 
as (To,Ko) = (1.35692,3.59389 x 10 13 ) where Ko is in cgs units, and the dividing density po used to model the transition 
from the crust to the nuclear matter is tabulated in the log of po [g/cm 3 ]. In the following calculations for BNS, a spherical 
solution of each EOS with gravitational mass Mi = 1.35M© is used as a reference, whose rest mass Md in solar mass units, 
circumferential radius R in km, compactness Mi/R in the geometric unit G = c = 1, and log of the central density p c in g/cm 3 
are tabulated. 
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FIG. 2: Selected components of hij along the a:-axis of the WL solutions with the orbital radius d/Ro = 1.5 for parametrized 
EOS 2H, HB, and 2B, from top to bottom panels of both sides respectively. Left (right) panels are Log-Linear (Log-Log) plots, 
where the a;— axis is normalized by A := n/Q. Upper and lower thin solid lines in the right panels are, respectively, the h yy and 



h zz components of the asymptotic solutions (|52[) of two point mass, Mi = 1.35M0 each, separated as the numerical solutions, 
d/Ro = 1.5. h+ in left panels is defined by h+ := (h xx — h yy )/2. 



up to the boundary of the computational domain. Up- 
per and lower thin black lines in the right panels are, 
respectively, the h yy and h zz components of the asymp- 
totic solutions (|52|) of two point masses. These lines do 
not exactly match the countours of the correspond- 



ing of numerical solutions for several reasons, including 
finite-size and higher order post-Newtonian effects. How- 
ever, the lines shift systematicaly from the numerical hij, 
which suggest that the numerical hij scales properly in 
the asymptotic region (as well as in the near zone) as the 
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compactness increases. 



B. Quasi-equilibrium sequences with different 
compactness 

A constant-rest-mass sequence of quasi-equilibrium so- 
lutions for irrotational BNS is considered as a model for 
the last several orbits of inspiral before merger. Such se- 
quences are computed for the models with different EOS 
parameters listed in Table 11111 The fixed rest mass of 
each model is that of a spherical star whose gravitational 
mass is Mi = 1.35M . Quantities of the spherical star 
for each model are also presented in the same Table. 

In Figf51 the binding energy E h := Madm — M and 
the total angular momentum J, normalized by twice the 
gravitational mass of the spherical star M — 2Mj , are 
plotted for models 2H, HB, and 2B. In the top and mid- 
dle panels, the results of the WL sequences are compared 
with the results of IWM sequences and of non-spinning 
point particles in 3PN circular orbits. Clearly, the IWM 
sequences coincide with the 3PN curve up to smaller sep- 
aration (larger f2M), whereas the WL sequences signifi- 
cantly deviate from the 3PN sequence. As the compact- 
ness (in this case from 2H to 2B) increases, the curves of 
the IWM sequence around the smallest separation come 
closer to the 3PN curve. In contrast to this, deviations 
of the WL sequences from the 3PN curve are even larger 
for the larger compactness. 

In the bottom panel of the Fig|3l the binding energy 
Eh := Madm — M of the WL sequences are compared 
with the results of the NHS sequences. Clearly, the dif- 
ference in the binding energy of two formulations is less 
than a percent; that is, the WL solutions almost coincide 
with the helically symmetric solution in the near zone. 

In [13, we have derived asymptotic conditions for 
equality Madm = Mk of the ADM and Komar masses 
[401 ] , which is related to the relativistic virial relation for 
the equilibrium [4lj . 



-0.005 



(54) 



In the WL/NHS formulation, the asymptotic fall-off of 
each field is sufficiently fast to enforce the equality. In 
Fig. [4j we evaluate the values of the fractional differences 
I Madm — Mr | /Madm for the WL sequences with the 
parametrized EOS 2H, HB, and 2B. The plots show that 
the differences are less than 2 x 10 -4 . The compactness 
increases in the order of 2H, HB, and 2B; the fractional 
differences, however, do not necessarily increase with in- 
creasing compactness in this range M\/R < 0.2. The fact 
that the fractional difference is well controlled for these 
sequences is evidence that the binding energy in Figj3] is 
calculated accurately. The virial relation Eq. (|54"|) , nor- 
malized by Madm, is also calculated to examine the ac- 
curacy of the numerical solutions, whose absolute value 
is about 0.5 ~ 1 times that of the fractional difference of 
two masses. 
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FIG. 3: Plots of the WL, NHS, and IWM sequences for the 
parametrized EOS 2H, HB, and 2B. Top panel: Binding en- 
ergy _Eb = Madm — M normalized by M = 2Mi with respect 
to the normalized angular velocity f2M of the WL and IWM 
sequences. Middle panel: Total angular momentum J normal- 
ized by M 2 of the WL and IWM sequences. Bottom panel: 
Normalized binding energy of the WL and NHS sequences. In 
each panel, a thin solid curve corresponds to that of the 3PN 
approximation. 
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FIG. 4: Fractional differences of Madm and Mk with re- 
spect to the orbital radius d/Ro of the WL sequences for 
parametrized EOS 2H, HB and 2B. 
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FIG. 5: Binding energy of WL sequences for the four 
parameter fitted EOS. Top panel: for the EOS, BGN1H1, 
SLy, and FPS. Bottom panel: for the EOS, ALF3, and APR1. 



C. Quasi-equilibrium sequences with 
four-parameter fitted EOS 

In the paper [l8j], optimal values for the parameters 
of four-parameter fitted EOS have been derived for 34 
candidates of the neutron-star EOS (17 selected EOS of 
nuclear matter with varied parameters). We choose five 
representative EOS, which are SLy [H, APR1 [H, FPS 
0, BGN1H1 [H, and ALF3 0. The first three are 
made only from normal nuclear matter, while BGN1H1 
involves a mixed phase with hyperons, and ALF3 with 
quarks. For the latter two EOS, the value of T becomes 
smaller in the mixed phase with the exotic matter at a 
few times above nuclear density [Hj]. However, BGN1H1 
is a stiff EOS having the largest p\ among them, and 
hence the core of the mixed phase is not large for the 
mass Mi = 1.35M . 

In Fig|5J the binding energy E\, of the WL sequences 
for these parametrized EOS are plotted. As in the case of 
the one-parameter parametrized EOS in Sec lIV Bl the se- 
quences with higher compactness Mi /R extend to higher- 
values of QM. Also, the WL sequences deviate from the 
3PN curve at larger QM. Among these EOS, APR1 is 
the softest, giving the most compact neutron-star model; 
and the corresponding binary sequence reaches the high- 
est value, ~ 0.058, of SIM. However, as seen in the bot- 
tom panel of Fig. the binding energy curve of APR1 is 
slightly off from the 3PN curve even for the smaller SIM 
of the sequence. In our neutron-star code, using a finite 
difference scheme, the core of the neutron star is cov- 
ered by fewer grid points in the central coordinates when 
the binary separation becomes larger and the neutron 
stars more compact; this may increase the numerical er- 
rors. We plan to incorporate a binary computation in the 
new code [37| . in which enough grids are maintained, to 
densely to cover the neutron star, irrespective of the bi- 
nary separation or neutron-star radius. The results of the 
APR1 curve as well as more compact binary sequences 
will be studied using the new code. 

In [l5| , the gravitational waveform computed from in- 
spiral simulations has been analyzed to estimate the ac- 
curacy with which gravitational wave observations can 
constrain neutron-star radius, an EOS parameter corre- 
lated with the departure from point-particle inspiral. A 
promising result is that the neutron-star radius can be 
constrained to SR ~ 1km for an interferometric detector 
with the sensitivity of Advanced LIGO, in either a broad- 
band configuration or a narrowband with peak sensitivity 
around 1150Hz. This suggests that the successful obser- 
vations of gravitational waves may exclude even a couple 
of EOS shown in Figf5] 



D. Comparison of the orbital phase in the last 
several orbits 

In this section, we approximately determine the orbital 
evolution in the late inspiral phase up to the onset of 
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merger using the quasi-equilibrium sequences computed 
in the previous section. To construct a quasi-equilibrium 
sequence, one assumes that each BNS evolves adiabati- 
cally along the sequence, that the radial velocity is much 
smaller than the orbital velocity. Given the rest mass and 
the EOS, each quasi-equilibrium sequence is defined by 
one parameter: The total energy and angular momentum 
of the binary system along a sequence are parametrized 
by the orbital angular velocity as E(fl) and J (ft). 

The time evolution of the angular velocity then be- 
comes 



dt \dn dt v ' 



(55) 



For the gravitational wave luminosity, dEjdt, we adopt 
the 3.5PN formula for two point masses [47]] ■ Tidal de- 
formation of the neutron stars in close orbits makes the 
attractive force between two stars stronger, and hence 
it accelerates the orbital velocity, resulting in the en- 
hancement of the gravitational wave luminosity. Thus, 
the 3.5PN formula for the luminosity is likely to under- 
estimate that of the BNS. However, this effect plays an 
important role only for the last ~ 1 orbit, and for most 
of the late inspiral orbits, the 3.5PN formula is a good 
approximation. 

Numerical integration of Eq. (|55p provides the relation 
between t and from 



t = J dQF(Q). 



(56) 



From this, the angular velocity as a function of time, 
fi(t), is obtained. Using this relation, we can also com- 
pute the approximate orbital phase evolution by 



N 



1 

2^ 



n(t)dt. 



(57) 



We note that the numerical model with the maximum 
value of f2 for each sequence presented in this paper does 
not exactly, but does approximately, correspond to a so- 
lution at the closest orbit. We stop the integration of Eq. 
(|56|) when f2 reaches its maximum. 

In the top panel of Fig. [5J VlM is plotted as a func- 
tion of time for EOS 2B, HB, FPS, and SLy in the WL 
formulation. In the bottom panel of Fig. [6j the results 
for 2B and HB, calculated in both the WL and IWM for- 
mulation, are compared. We also plot the results of two 
point masses, derived from the Taylor-T4 formula 48]. 

The top panel of Fig. [S] shows that for the small val- 
ues of f2, all the curves approximately agree, irrespective 
of the EOS. This is natural because for such small val- 
ues, tidal deformation does not play an important role 
and orbital velocity is sufficiently small (v < 0.3c) that 
the post-Newtonian formula (Taylor-T4 formula) with 
the point-particle approximation should be an excellent 
approximation. 

By contrast, the values of £l(t) computed from the nu- 
merical sequences deviate from those given by the Taylor- 
T4 formula for fiM > 0.035-0.04, for all of the EOS 
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FIG. 6: Top panel: Orbital angular velocity, f2, as a function 
of time for EOS 2B, HB, FPS, and SLy of the WL sequences. 
Bottom panel: Same as the top panel but for EOS 2B and 
HB of the WL and IWM sequences. For the both panels, the 
results by the Taylor-T4 formula are also plotted. The units 
of fi and time are M~ and M, respectively. For comparing 
the results, the time axis is shifted such that QM = 0.03 is 
aligned at t = for all the curves. 



and all formulations used to compute the quasiequilib- 
ria. This is due to the tidal deformation of the neutron 
stars; the rate of change of the energy as a function of Q 
approaches zero for the close orbits, as seen in Figs. [3] and 
[5] This deviation occurs at more distant orbits for less 
compact neutron stars (i.e., for the stiffer EOS), indicat- 
ing, as expected, that one can extract from the curve f2(t) 
a characteristic of the component neutron stars related 
to their compactness and a corresponding parameter of 
the EOS. 

The bottom panel of Fig.[H]shows that the curves f2(i) 
computed by the WL and IWM formulations are signif- 
icantly different, as expected from the results of E(il). 
In the case that the IWM formulation is adopted, the 
merger time is overestimated by ~ 50M, which is a quite 
a large factor. This suggests that the results in the IWM 
formulation do not work well for predicting the evolution 
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FIG. 7: Top panel: Orbital cycle, N, as a function of time for 
EOS 2B, HB, FPS, and SLy of the WL sequences. Bottom 
panel: The same as the above but for EOS 2B and HB of the 
WL and IWM sequences. For the both panels, the results of 
the Taylor-T4 formula are also plotted. For comparing the 
results, the time axis is shifted such that N = is aligned at 
flM = 0.03 for all the curves. 



of the last several orbits before the onset of merger. 

In Fig. [71 we plot the curves of N as a function of QM; 
the top panel is for EOS 2B, HB, FPS, and SLy in the WL 
formulation and the bottom for EOS 2B and HB in the 
WL and IWM formulation. The top panel shows that the 
number of orbital cycles in the late inspiral phase depends 
strongly on the EOS. For a soft EOS, e.g., EOS 2B, in 
which the compactness of the neutron star is largest, the 
number of cycle is largest. By contrast, for a stiff EOS 
such as SLy, the number of cycles may be smaller by ~ 1 
than that for EOS 2B. 

In bottom panel of Fig. [7J the results for the number 
of cycle calculated from different formulations are com- 
pared. As expected from the results for fi(t), the IWM 
formulation overestimates the number of cycles. The er- 
ror AN is ~ 0.5 for the EOS 2B; i.e., one cycle of gravi- 
tational waves would be overestimated. 



V. DISCUSSION 

The deviations of the binding energy and total angular 
momentum of WL/NHS sequences from the 3PN point- 
particle sequence as well as from the IWM sequences are 
likely to be due to the tidal deformation of neutron stars 
in the binary system coupled with general relativistic 
effects. As the compactness of the component neutron 
stars increases, the deviation from the 3PN sequence at 
a certain value of fiM decreases - WL/NHS sequences 
become closer to the point-particle sequence, but not by 
as much as the IWM sequences do. It has been believed 
that, as the compactness of the component neutron stars 
increases, the behavior of the binding energy and angu- 
lar momentum of binary sequences more closely approx- 
imates that of point masses. This is found in the results 
of IWM sequences but to a lesser extent in the WL/NHS 
sequences. The behavior of the IWM sequence was in- 
terpreted as the effacing of the tidal effects due to the 
strong gravity: that is, as the compactness increases, the 
sequences of binary neutron-star solutions become much 
closer to the sequences of two point masses, because the 
tidal effect is masked by the stronger self-gravity of each 
component star. However, the results of WL/NHS se- 
quences suggest that such effacing of the tidal effect seen 
in IWM sequence is an artifact of the conformally flat 
approximation, at least for the case of equal mass binary 
neutron stars. 

In the WL/NHS formulations, all components of Ein- 
stein's equation are solved without approximation on a 
initial hypersurface, while in the IWM formulation, some 
terms of second post-Newtonian order are truncated. As 
discussed in Q the difference between the IWM and 
WL/NHS formulations in the binding energy Ej, is es- 
timated at second post-Newtonian order as Mh a bV a v , 
where the magnitude of the orbital velocity v a is typ- 
ically v ~ 0.34(f7M/0.04) 1 / 3 . Since h ab is 0(v 4 ), the 
order of the difference in the binding energy is given by 
AEb/M — 0(v 6 ) ~ 10 -3 , and a larger deviation as v 
becomes larger for more compact sequences is expected. 
This estimate is consistent with our results shown in 
Fig. [3] and O Note also that the tidal effect is larger for 
the EOS with a larger T as we used in our computations. 
So far, our WL/NHS codes have passed several code tests 
(as have the IWM codes), and results of two independent 
WL codes agreed for a BNS sequence with M\JR = 0.17 
as shown in [7j . These results support our argument that 
the WL/NHS results accurately correct the IWM results. 
A computation of quasi-equilibrium BNS sequences us- 
ing a totally different numerical method, such as the fully 
constrained scheme [9], would be a helpful additional 
check. 

We think our results suggest that the circularity of 
orbit is more accurately enforced on a WL/NHS se- 
quence than a IWM sequence. However, in such quasi- 
cquilibrium sequences, some important features of real- 
istic inspirals are ignored. Those include the radial ve- 
locity due to gravitational radiation reaction at 2.5PN 
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order that is likely to be comparable to the 2PN terms 
during the last few orbits where the neutron-star velocity 
is of order v ~ 0.1, and a tidal lag angle of about 10 — 20 
degrees that is found in inspiral simulations. Therefore, 
a caveat is that estimates of the merger time and or- 
bital cycles using quasi-equilibrium sequences shown in 
Sec. IIV Dl involve errors due to ignoring these effects. 

Recently, several groups have developed methods to 
treat the general relativistic tidal deformations analyti- 
cally [49|]. Comparison of these analytic results and the 
present results for WL/NHS sequences may be useful in 
calibrating the binding energy or the total angular mo- 
mentum of the quasi-equilibrium sequence in the regime 
where the relativistic tidal effects become important. Fi- 
nally, by combining the analytic and numerical results, 
more accurate quasi-equilibrium models for the late in- 
spirals may be constructed [50| . 

The WL/NHS formulations can be also used to con- 
struct models of rotating neutron stars. In [5l|, ax- 
isymmetric rotating relativistic stars are computed using 
the fully constrained formulation with maximal slicing 
and the generalized Dirac gauge conditions Q. Those 
solutions agreed with the ones calculated using a sta- 
tionary axisymmetric metric with the additional discrete 
symmetry of the simultaneous transformation, t — ► — t 
and <j> — ► — (j). The WL/NHS formulations include more 
general stationary axisymmetric spacetimes, which do 
not depend on the additional symmetry. Therefore, the 
WL/NHS formulations can be applied, for example, to 
rotating neutron stars that may have both toroidal and 
poloidal components of the magnetic fields as well as 
meridional circulation. Even in this case, the WL/NHS 
formulation can be used to compute exact equilibria that 
are more general than those calculated in [52|]. We plan 
to extend our codes to compute relativistic rotating stars 
and binary systems that each include strong magnetic 
fields. 
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APPENDIX A: BASIC EQUATIONS 

In this Appendix, the system of equations used in 
WL/NHS codes is presented in detail. The equations 
include all components of Einstein's equation, the first 



integral of the relativistic Euler equation, and the rest 
mass conservation equation for the irrotational flow. The 
WL/NHS formulations are based on 0, El] • 



1. Conventions 

As mentioned in Sec lII Al the 3+1 decomposition is ap- 
plied to the spacetime M. in the WL/NHS formulations. 
First, several definitions for the quantities relating to the 
spatial geometry are introduced. 



a. Connections 

The spatial metric j ab , a conformally rescaled spatial 
metric j a b, and a flat metric f ab are associated with the 

o 

derivatives D a , D a , and D a , respectively. We introduce 
the conformal rescaling by j ab = ip 4x f a b, whose determi- 
nant 7 is equal to that of the flat metric /, 7 = /, to 
specify the decomposition of the spatial metric uniquely. 
Covariant derivatives D a and D a are related by 

D b X a = D b X a + CS C X C , (Al) 

where X a is a spatial vector, and a coefficient C% b is 
written 

Clb = \l cd {Daldb + D blad -b dlab ) 

= ^{l C bDa^ + l C aDbi>-labl Cd D d ^). (A2) 
o 

Also, D a and D a are related by 

D b X a = D b X a + C b a c A c , (A3) 
where C^ b is written 

C c ab = ^l cd (Daldb+D b j ad - D d % b ) 

= \l cd {D a h db + D b h ad -D d h ab ). (A4) 

A trace of C° b 

C b ba = \l bc D a j bc = ^4V7, (A5) 

and the condition 7 = / that specifies the conformal 

~ o 

decomposition imply C^ a = and hence D a tjj — D a ip. 
The relations 

~ab c c b + -bc^ + f, a ~ac = Q (Ag) 

o 

and 7 = /, and the Dirac gauge condition D b Y = 
imply r h C r ab = 0. 
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b. Conformally rescaled extrinsic curvatures 



c. Conformally rescaled intrinsic quantities 



The form of the extrinsic curvature K a b is discussed in 
Sec lII Bl In the equations for our numerical code, it is 
decomposed in terms of the trace K = j ab K a b and the 
tracefree part A a b, 



Kab = A ab + -j a bK. 



(A7) 



The conformally rescaled tracefree part A a b is defined as 

A b = A a \ (A8) 

and its index is lowered (raised) by ^ a b (7° b )- 

We define Lx^ab as the tracefree part of £xjab, where 
X a is a spatial vector on St, 



Lx7ab = £x7ab ~ 3 labl^ £^Xlcd 



(A9) 



= D a X b + D b X a - - lab D c X c (A10) 

The r.h.s. of this equation is a conformal Killing operator, 
and its conformally rescaled version is defined by 



Lxlab = tp^Lxlab- 



(All) 



Note that a vector is rescaled, X a = X a , and ^ a b is used 
when lowering the index of the rescaled vector. 

When helical symmetry, £kg a /3 — 0, is imposed as in 
Sec lII Bl the tensors A a b and A a b have the forms 

Aab = -^-L^jab and A ab = — A/y &, (A12) 
la la 

respectively; while for the WL formulation, 

Aab = -^-Lpjab and A a b = ^-Lplab- (A13) 
la la 

The following expression for the conformally rescaled A a b 
is used later, 



1 



A a ° = ^ [D a P» + D b p a ~ ^la h D c F 



(A14) 



The Ricci tensor 3 R a b of the spacelike hypersurface St 
associated with the spatial metric j a b is decomposed into 
terms related to the conformal factor ip, 3 R^ b , and the 
conformal Ricci tensor 3 R a b associated with jab- 



3p _ 3r>V , 3n 
n a b — ti a b + K ab- 

The first term is written 



(A17) 



3 Rf b = -^D a lhr- U-WD.-r 



2 - ~ , „ 2 

+ ^b aV jD bV j-^ ab ^b c ^b c ^. (Ais) 

In 3 R a b, terms linear in h a b or h ab are separated as 

3 R ab = - \b C I) c hab + Rab + Rat, (A19) 



where R^ b includes terms linear in the conformal metric 



in the form of flat divergences 



R^b = -\(faAF c + hcD a F c ), (A20) 

~ab _ n_ uab. (A21) 



F a := DbT = D b h 
non-linear terms, R^ b h , are written 



Rat = -\{D b h cd D c h ad +D a h cd bMd + h cd D c bdhab) 

yk r^c 1 r^c /~id r^d (~ic 
U a^ c b T ^ a b^dc °ac U bd 

-\[b b (.h ac F c ) + D a (h bc F c ) ] + F c C cM , (A22) 

where C c ^ a b ■— lcdP d b - The above expression for R^ b 
can be simplified by applying the condition 7 = / and 
the generalized Dirac gauge condition, implying C ba = 
and F a = 0. 

The Ricci scalar curvature 3 R of S< is related to the 
conformal Ricci scalar 3 R := ^ ab3 R ab by 



(A23) 



The last term in the above 



L^Jab = £(/>7ab — T^lab^ £<plcd 



(A15) 



= D a & + D b <f )a --%bD c <f ) c , (A16) 

with cf> a = (ff and cj) a := ; y a b l fi b , appears only in the he- 
lically symmetric case and is eliminated when the WL 
formulation is used. 



2. Equations for the gravitational fields 

Equations used in the numerical code are shown below. 
Although we impose the gauge conditions (fl"2"|) and (|13p , 
the following equations are not restricted to these choices. 
The conformal decomposition, however, is specified by a 
condition 7 = / that is used, for example, to obtain the 

~ o 

relation D a %p = D a ip. 
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a. Hamiltonian constraint 

The projection of Einstein's equation along the normal 
n a to the hypersurface yields 

{G af} - %-KT aP )n a nP 

= U 3 R + K 2 - K ab K ab - 16ttph) = 0. (A24) 



Substituting Eq. (|A23[) . we have 
(G af3 - 8nT a0 )n a nP = -i 



0. (A25) 



The above equation is rewritten to isolate the flat 

O 

Laplacian A^ := D a D a %p on the l.h.s., and the other 
terms are treated as a source on the r.h.s., 



with the source Sn given by 



(A26) 



S H = -h ab D a D b iP + j ab C° b D c ip + 

o 

( A ab A ab - 2 -K 2 \ - 2n^ Pu . (A27) 



b. Momentum constraint 

The momentum constraint is written in an elliptic 
equation to be solved for the covariant component of the 
conformally rescaled non-rotating shift (3 a := j a bP b - We 
begin with 



{G aP - 8^T Q/3 ) 7a a r/ 

= -D b K a b + D a K + 8irj a 



1 - 



2 ~ 



~D b [ip 6 A a b j + ^D a K + 87rj a = 0,(A28) 



then substitute Eq. (IAl4j) and a relation 

D b D a j3 b - D a D b p b = 3 R ab p b , (A29) 

to obtain 

D b D b f3 a + ^D a D b (3 b + 3 R ab f3 b + nb b (L^) a b 



+2aA b — D b (^J - - aD a K - 16naj a = 0. 



(A30) 

From the first two terms of the r.h.s. of Eq. (IA30|) . the 

o - loo- 

flat terms A(3 a + -D a D f3 b are similarly isolated, 
D b D h j3a + \DaD b p b = kj3 a + ^D a D b p b 



OO- 



+h bc D b D c (3 a - j bc D b (CiM - j bc CtDdPa 
~f c CtDc~Pd + \D a (h bc D b p c - f c Ctp d ). (A31) 



We keep D a instead of replacing it by D a and a connec- 
tion C c ab in a couple of terms in the Eq. (|A31[) , to shorten 
the equation. A decomposition proposed by Shibata, 

fl, = G a + \i) a (B - x b G b ), where D a x b = S b a , (A32) 



is substituted in the expression for the flat operator 
A/3 a + ^D a D b p b , 

U a + \£) a D b p b = AG a + ^-D a ( AB — x b AG b ), (A33) 
3 6 

to obtain elliptic equations that are solved simultane- 
ously, 



AG a = S a , 
AB = x a S a , 

where the source S a is written 



(A34) 
(A35) 



S a := -h bc D b Dj a + j bc D b (C^Jd)+l bc Ci c D d p a 
+f c CtDj d - ^D a (h bc D b p c ~ f c Ctj d ) 



3 R ab p b - nb b L^ ab - 2aA a °^D b 



a - 



-- aD a K + Wiraja- 



(A36) 



A term D b L,pj ab is computed from 



D b L^ ab = t c D c L,t,% b -C^ bd L^ cd 

+ Dcl cb L^lab + C d dc T b hlab, (A37) 
which is dropped when the WL formulation is used (see, 

secnnji. 

c. Spatial trace part of Einstein's equation 

The spatial trace of Einstein's equation is combined 
with the Hamiltonian constraint, 

(G Q/3 -87rT a/3 )( 7 ^ + Vr/) 

= - - 3 R + -D a D a a + 2£ n K 
4 a 

-\{K 2 + 7K ab K ab ) - 4n(p H + 25) = 0, (A38) 

and it is solved for the combination atp. Using a relation, 

lo 2 
- - 3 R + -D a D a a 
4 



o 



aip 5 



D a b a {a^) -^ 3 R 



(A39) 
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and applying helical symmetry, the above equation is 
rewritten 

D a D a {aiP) - ^ 3 R - tP 5 £ u K 



aip 5 



nr r - (lA ab A ab + ^K 2 ) - 2r,-n, •'(,,„ +2,S'l 



= 0. 



(A40) 



Isolating the flat part A(at/j), an elliptic equation is de- 
rived 



A(e*V>) = S tI , 
where the source <St r is written 



(A41) 



Str := - h ab D a D b (aij) + r b C c ab D c (a^) + ^ 3 R 

o 

+ 2 7 rai/; 5 ( / 9H + 25). (A42) 



d. Spatial tracefree part of Einstein's equation 



The tracefree part of the tensors are also denoted by sub- 
scripts TF, hereafter. We further eliminate terms propor- 
tional to 7afc remaining in this expression for £ ab later in 
this section. 

We derive two different equations to solve for h ab . One 

o 

is an elliptic equation in which Ah ab is separated from 
3 R ab as in Eq. (|A"l9)> ; it is used for both the WL/NHS 
formulation. The other is for the NHS formulation in 

o 

which an operator (A — Q, 2 d^)h ab is separated. The <f> 
derivative term in this operator is separated from a term 
- a -£ u >K ab , which is derived by applying helical symmetry, 
(HHJ) and (fl"T)) . to the time derivatives. For the former 
equation, the above £ a b is rewritten 

£ ab = ~Ah ab + RV b + Rl<t + 3 Ri b --D a D b a 
i a 

+KK ab - 2K ac K b c + -£^K ab - 8^5 afc , (A47) 
a 

and for the latter, 



£ ab = -\ (A - n 2 d 2 ) h ab + r» + kl 



NL 
6 



- 3 Rt b - -D a D b a + KK ab - 2K ac K b c 
a 

— £uK ab — -Q 2 d\h ab — 8nS ab . 
a 2 v 



(A48) 



The projection of Einstein's equation to the initial hy- 
persurface E t is written 

= -£ n K ab + j ab £ n K + 3 R ab — - jab 3 R 

+KK ab - 2K ac K b c - i 7afc (K 2 + K cd K cd ) 

— (D a D b a-j ab D c D c a)-8nS ab . (A43) 
a 

The equation to solve for the non-conformal part of the 
spatial metric h ab is derived from the tracefree part of the 
above equation (|A43j) ■ The tracefree operation eliminates 
terms proportional to j ab . Applying helical symmetry, 
(fir]]) and (IT]), the tracefree part of Eq. (|A43|) is written 



{Gap - 87rT a/S )( 7o <V - - lab T P ) = £ T J = 0, (A44) 

where £ ab is defined by 

1 3 1 

£ a b '■= —£ujK ab + R ab D a D b a 

a a 

+KK ab - 2K ac K b c ~ 8irS ab , (A45) 
and £j b F is its trace free part 



pTF _ 

c ab ■— 



(ja C lb d -^JabJ Cd )£cd = (ja C Jb d -^JabY d )£cd- 

(A46) 



Terms proportional to j ab in Eqs. (|A47p and (|A48p are 
now eliminated further to simplify the equations. Intro- 
ducing barred quantities, 



2 - - 



6 ^ 



-D a D b iP + —}D a ipD b il), (A49) 



D a D b a = D a D b a - C c ab D c a 



--(D a aD b iP + D b aD a ^), (A50) 



their combination becomes 



3 Rt h - -D a D b a 



-l^D a D b (a^ 2 ) + J-C c ab D c (a^ 2 ) 
aip* dip 

D a {aip)D b ip + D b (on/j)D a ip\ , (A51) 



aip 



which satisfies 



/ _ , | \ TF f _ . 1 \ TF 

( 3 Rtb--D a D b a) = [ 3 Rt b --D a D b a) . (A52) 

Next, substituting K ab = A ab + ^ a bK to terms relating 
to K ab , their tracefree part satisfies 



KK ab - 2K ac K b c + -£^K ab ] TF 
a ) 
1 1 \ tf 

-KA ab -2A ac A b c + -£ UJ A ab ) . (A53) 
A a J 
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For the matter source term, 

Sab = TafiJ^lb = (e + p)u a u b + j abP , (A54) 
where u a := "f a a u a . We also introduce a barred quantity 
Sab ■= phu a u , (A55) 
that satisfies Sj b F = Sj b F , where h = (e+p)/p is used. 

o 

The tracefree operation to the operator (A — Sl 2 <9| )h cc i 
is written 



\%b^) ( 



a - n z di h cd 



a - n 2 dl h a 



+ \labD e h cd D e h cd - ^% b n 2 d^h c %h cd 



,(A56) 



where relations j cd D e h cd = 7 cd d^h cd — implied by 
7 = / is used. The same operation to the Laplacian is 
written similarly as above, but without 8$ terms. 

Finally, the trace free part £j F = results in the fol- 
lowing elliptic equation, 



Ah a b — S a b 

where the source S a b is defined by 



1 



Sab := 24V - -i ab D e h cd D e h cd) (A58) 

and £j b F is a tracefree part of £ ab , which is written using 
the rescaled A ab , 



Sab := Rab + Rab + ^ab ~ ~D a D b a 

1 



SNL I 361/' 



1 



+ ^ 4 KA ab - 2ip 4 A ac A b c 

+ -£ LU (^ 4 A ab )-8nS ab . (A59) 
a 

o 

For the equation with the operator A— fi c?3, it is written 
(k-tfdfjKb^Sab (A60) 

with 

Sab ~ 2£~Ib F ~ \labD e h cd D e h cd + ~% b n%h C %h cd . 



(A61) 



Using the rescaled A a b, S a b is defined by 

1 



c dD i pNL , 3p0 

Cab ■— n ab -+- n ab -f n ab 



1 



-D a D b a 



+ -^KA ab - 2^ 4 A ac A b c 

+ -£^ 4 A ab ) ~ \tfdlhab - 8TTS ab , (A62) 

where a difference from (|A59I) is a term in the last line. 



e. Matter source terms 

In the above, the matter source terms, pn, j a , S and 
S ab , that appear in the field equations are obtained from 
the stress energy tensor. We write the projection of the 
stress energy tensor in terms of the fluid variables and 
metric potentials. The 4- velocity for irrotational flow 
u a = u*(fc" + v a ) is decomposed with respect to the 
foliation St as 

u a n a = -au l (A63) 

u a laa = u a = \D a ^> = \f) a <S>, (A64) 
h h 

where the velocity potential $ is introduced by hu a — 

Using these relations, the matter source terms of the 
field equations become 



p H := T a/3 n a n = hpiau 1 ) 2 ~ p, (A65) 

(A66) 



j a ■= -T aii "ia a n (i = pavfDa®, 



S := T a ^ aSi = hp[(c 



l]+3p, (A67) 



Sab := T al3l a a jb = J^Da^Db® + Plab, (A68) 



(A57) 

or with a barred quantity, 



S ab = $D a $D b $. 
h 



3. Equations for irrotational fluid 



(A69) 



Following Sec lIICI and lllDl a set of equations used in 
our codes to solve for the matter variables are derived. As 
independent variables, we choose the relativistic enthalpy 
per baryon mass, the time component of the 4- velocity, 
and the velocity potential, {ft,, it', $}. For the first two 
variables, the first integral Eq. (|34[) and the normalization 
of the 4- velocity u a u a — — 1 are solved. Using a relation 
derived from Eqs. ([3T|) and 



(A70) 

, (A71) 
(A72) 



v a +uj a = -i-£) a $, 
h/ir 

these equations are rewritten, 
1 



a 
1 

a 2 h 



(£ + u a D a <Z>y - D a <Z>D a <f> 



1/2 



— (£+UJ a Da<S>) 



T, and the second 



where the first one is from u a u a 
from Eq. (|3"4]) . 

An equation to calculate the velocity potential <£> is 
derived from the rest mass conservation law, Eq. (|2^|) . 



1 



-.£ u (py/^g) 



1 



■.^(pU^y/^f) 



-Daiap^v 11 ) = 0. (A73) 
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Substituting Eq. (|A70|) in the above relation, we have an 
elliptic equation for <!>, 

D a D a $ = D a (fcuVV(A»$-feu t w«)— £°^r- (A74) 

ap h 

This equation is solved with Neumann boundary condi- 
tion to impose the fluid 4- velocity u a to follow the surface 
of the star. The surface is defined by the vanishing pres- 
sure p = 0, which coincide with the h — 1 surface in our 
EOS (see, Sec III Dp . Hence, the boundary condition is 
written 



u a V a h = at h=l. 



(A75) 



and, using £kh = and Eq. (|A70|) . Neumann boundary 
condition for the potential $ is rewritten, 



{D a <S> - hu l uj a )D a h = 0. 



(A76) 



where D a h is normal to the stellar surface. 

Finally we rewrite the above set of equations for the 
helically symmetric irrotational flow using the flat deriva- 

o 

tive D a ; 



h = 



1 / ° \ 2 1 

( £ i r.a n *l z,ab i 

a 2 V 



£ + u> a D a $) -—j ab D a $D b $ 



i> 4 



1/2 



A$ = S, 



(A77) 
(A78) 
(A79) 



where S is defined by 



S = - h ab D a D b <t> + r h C c ab D c ® - ^r b D a iPD b <S> 
+ -^Cj a D a (hu t ^, 6 ) + ^ i hu t D a u a 



ap 

ap h 

and, for 7 = /, D a u a = D a Cj a = D a f3 a . 



(A80) 



component satisfies a field equation whose principal part 
is the same as that of a scalar equation. 

For the spatial tracefree part of Einstein's equation 
solved for the non-conformally flat part h ab , writing the 

o o 

principal part C :— A or A — J7 2 9|, the field equations 
become 



£h ab — S ab 



(Bl) 



Expanding each cartesian component of h ab in scalar 

o 

multipoles, the equation with the operator A — fi 2 3| 
becomes a Hclmholtz equation for each mode, 



A + mW)hXY em =S t a TYt m . 



(B2) 



Hence these elliptic equations are integrated using 
Green's formula, 



h ab (x) =--?-/ G(x,x')S ab (x')d 3 x' 

47T Jy 



1 

47T 



dV 



G{x,x')D ,c h ab (x') 
-h ab {x')D' c G(x,x')\ dS' c . (B3) 



where x and x' are positions, x,x' € V C E t , and the 
Green function G(x,x') satisfies 



CG(x,x') = -An5(x-x'). 



(B4) 



We choose the Green function G(x, x') without boundary 
for the BNS calculations. 

o 

For the Laplace operator, C = A, a multipole expan- 
sion of G(x,x') in associated Legendre functions on the 
spherical coordinate is written 



G(x,x') 



1 



{l-m)\ 
(£ + m)\ 



1 1 e=o m=0 

xP e m {cosd) P i m (canO f ) cos m(cp - ip'), (B5) 

where the radial Green function gt(r,r') becomes 



APPENDIX B: SELF-CONSISTENT FIELD 
ITERATION SCHEME 



gi(r,r') = 



(B6) 



1. Elliptic equation solver 



As mentioned in Sec lIIIBl components of the metric 
are computed on a spherical-coordinate grid whose ori- 
gin is placed at the center of mass. The momentum con- 
straints and the tracefree part of Einstein's equation are 
a spatial vector and a tensor equation, respectively, and 
it would be natural to write the equations in components 
along the spherical coordinates [36j . It is simpler, how- 
ever, to solve these equations for cartesian components, 
yet on the spherical coordinates, because each cartesian 



with r> := sup{r, r'}, r< := inf{r, r'}, and the coeffi- 
cients e m are equal to eo = 1 for m — 0, and e m = 2 for 
m > 1. 

o 

For the case with the Helmholtz operator, C = A + 
TO 2 f2 2 , we choose the Green function for the half-retarded 
+ half-advanced field p^j . 



G(x,x') = ^2 9e m (r,r')e r 



(£-m)\ 
(£ + m)\ 



1=0 m=0 

xP t m (cos 6) P t m (cos 9') cos m(tp-<p'), (B7) 
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where the radial Green function gi m {r, r') is constructed 
from the spherical Bessel function of the first and second 
kinds je(x) and ng(x), 



d/R 2d/M Ro/M QM Madm J/M 2 



l {r,r') = 



for m = 0, 



(B8) 

— mQ (2£ + 1) ji(mn r<) ni^mQ r>), 
for m > 1. 



2. Summary for iteration scheme 



Eq. (|B3[) is used as an elliptic equation solver for the 
field variables {a, /3 a ,tp, h a b}. In the code, the ellip- 
tic solver is used to compute the combination aip from 
Eq. (|A41|) ; to compute the potentials (|A32p of the shift 
vector f3 a from Eqs. (|A34j) and (|A35|k and to compute 
the gauge potentials (j4"9"j) from Eq. (|5T 



For the fluid variables, {h, u*,<I>} are found from 
Eqs. (|A77|) , (|A78|) . and (|A79|) . respect ively A detailed 
description of a method to solve Eq. (|A79p is found in 
[3^ |. As we use the surface- fitted coordinates to cal- 
culate neutron stars, the surface R(9f,4>f) becomes an 
additional variable. A stellar surface is defined by the 
pressure p = 0, and, instead, it is located by a condition 
q = p/p = in the code. 

A solution is specified by two parameters for an equal 
mass binary, which we take to be the orbital angular 
momentum and the injection energy, {£},£}. We in- 
troduce one more parameter Rq to normalize the ra- 
dial coordinate, where Rq is half the coordinate diam- 
eter of a neutron star along the (9f,<fif) = (tt/2,0) 
line. These parameters are calculated from the condi- 
tions R(ir/2, 0)/R = 1 and R(ir/2, n)/R = 1, after pre- 
scribing a value of a thermodynamic variable at a point 
in a star, for which a central value of h is fixed at r/ = 0. 
These conditions are applied to Eq. (|A77[) , and solved for 
the three parameters. 

All these variables are assigned on each grid point, and 
the parameters are calculated from the equations men- 
tioned above in each iteration cycle. If we represent the 
set of fluid and metric variables by 'J, we can describe 
the iteration schematically as follows. The variables are 
are updated from their values at the Nth iteration cycle, 
*( N) , to the (N+l)th, *( N+1 ), using softening, in the 
manner 



*( N +D = A* + (1-A)^ N \ 



(B9) 



where A is the softening parameter, chosen to be in the 
range 0.1 to 0.3 to accelerate convergence. For a criteria 
to determine convergence, a relative difference of succes- 
sive cycles 



2|^(N+l) _ ^(N)| 



|^(N)| 



< 6 



(B10) 



1.3125 
1.3438 
1.3750 
1.4375 
1.5000 
1.6250 
1.7500 
1.8750 
2.0000 
2.5000 
3.0000 



9.2784 
9.4190 
9.5654 
9.8715 
10.192 
10.873 
11.590 
12.330 
13.090 
16.234 
19.455 



7.0693 
7.0095 
6.9566 
6.8671 
6.7948 
6.6912 
6.6229 
6.5759 
6.5451 
6.4936 
6.4852 



0.030149 
0.029649 
0.029126 
0.028019 
0.026901 
0.024693 
0.022648 
0.020794 
0.019118 
0.014113 
0.010889 



2.67191 
2.67206 
2.67223 
2.67265 
2.67311 
2.67422 
2.67537 
2.67648 
2.67756 
2.68130 
2.68404 



0.96637 
0.96819 
0.97048 
0.97582 
0.98279 
0.99905 
1.0178 
1.0375 
1.0572 
1.1415 
1.2229 



TABLE IV: Solution sequence for the EOS 2H. 



d/Rg 2d/M Ro/M QM Ma 



DM 



J/M 2 



1.3750 
1.4062 
1.4375 
1.5000 
1.6250 
1.7500 
1.8750 
2.0000 
2.5000 
3.0000 



6.8380 
6.9403 
7.0511 
7.2800 
7.7600 
8.2681 
8.7975 
9.3425 
11.606 
13.931 



4.9731 
4.9353 
4.9051 
4.8533 
4.7754 
4.7246 
4.6920 
4.6712 
4.6423 
4.6435 



0.045877 
0.045086 
0.044231 
0.042526 
0.039177 
0.036025 
0.033145 
0.030526 
0.022621 
0.017521 



2.66507 
2.66521 
2.66536 
2.66580 
2.66688 
2.66810 
2.66932 
2.67057 
2.67511 
2.67860 



0.89082 
0.89216 
0.89378 
0.89812 
0.90829 
0.92072 
0.93485 
0.94926 
1.0132 
1.0789 



TABLE V: Solution sequence for the EOS HB. 



d/R 2d/M R /M MM Madm J/M 2 



1.4375 
1.4688 
1.5000 
1.5313 
1.5625 
1.6250 
1.7500 
1.8750 
2.0000 
2.5000 
3.0000 



5.5971 
5.6801 
5.7713 
5.8622 
5.9542 
6.1443 
6.5438 
6.9613 
7.3930 
9.1954 
11.050 



3.8936 
3.8673 
3.8475 
3.8284 
3.8107 
3.7811 
3.7393 
3.7127 
3.6965 
3.6782 
3.6835 



0.059912 
0.058871 
0.057759 
0.056652 
0.055575 
0.053414 
0.049241 
0.045416 
0.041907 
0.031170 
0.024219 



2.66109 
2.66118 
2.66133 
2.66150 
2.66171 
2.66211 
2.66318 
2.66432 
2.66555 
2.67036 
2.67433 



0.85642 
0.85733 
0.85840 
0.85959 
0.86099 
0.86424 
0.87225 
0.88246 
0.89311 
0.94283 
0.99755 



is used, with <5 = 10 6 in the present calculations. 



TABLE VI: Solution sequence for the EOS 2B. 



APPENDIX C: FORMULAS FOR MASS AND 
ANGULAR MOMENTUMS 

Definitions of the quantities shown in tables and fig- 
ures, which characterize a solution of BNS, and their 
expressions used in actual numerical computations, are 
summarized in this Appendix. 

The rest mass is the baryon mass density measured 
by comoving observers integrated over the initial hyper- 



20 



d/Ro 


2d/M R /M 




SIM 


A/adm 


J/M 2 


1 


.3750 


6.7360 


4. 


.8989 





.046804 


2 


.66479 


0.8804 


1 


.4375 


6.9455 


4. 


.8316 





.045130 


2 


.66505 


0.8987 


1 


.4687 


7.0572 


4. 


.8049 





.044256 


2 


.66526 


0.8983 


1 


.5000 


7.1692 


4. 


.7795 





.043403 


2 


.66546 


0.8998 


1 


.5625 


7.4012 


4. 


.7367 


0. 


.041683 


2 


.66600 


0.8952 


1 


.6250 


7.6424 


4. 


.7030 





.039990 


2 


.66655 


0.9093 


1 


.7500 


8.1425 


4. 


.6529 


0. 


.036778 


2 


.66777 


0.9109 


1 


.8750 


8.6640 


4. 


.6208 





.033842 


2 


.66899 


0.9397 


2 


.0000 


9.2007 


4 


.6004 





.031172 


2 


.67025 


0.9413 


2 


.5000 


11.431 


4. 


.5722 





.023105 


2 


.67481 


1.082 


3 


.0000 


13.721 


4. 


.5737 





.017898 


2 


.67834 


1.031 



TABLE VII: Solution sequence for the EOS SLy. 



d/Rp 2d/M Rq/M SIM Madm J/M 2 

1.6875 5.8207 3.4493 0.057394 2.66144 0.85585 

1.7500 6.0069 3.4325 0.055113 2.66195 0.85933 

1.8125 6.1950 3.4179 0.052972 2.66247 0.86339 

1.8750 6.3877 3.4068 0.050902 2.66302 0.86770 

1.9375 6.5845 3.3985 0.048909 2.66363 0.87233 

2.0000 6.7839 3.3920 0.047008 2.66423 0.87706 

2.5000 8.4399 3.3760 0.035057 2.66900 0.92198 

3.0000 10.148 3.3826 0.027269 2.67308 0.97231 

TABLE VIII: Solution sequence for the EOS APR1. 



d/Rp 2d/M Rq/M SIM M ADM J/M 2 



1 


.4375 


6.3451 


4. 


.4140 





.050851 


2 


.66331 


0.87458 


1 


.4688 


6.4477 


4. 


3899 





.049882 


2 


.66346 


0.87609 


1 


.5000 


6.5504 


4. 


3669 





.048919 


2 


.66364 


0.87782 


1 


.5313 


6.6534 


4. 


.3451 





.047977 


2 


.66387 


0.87955 


1 


.5625 


6.7597 


4. 


3262 





.047029 


2 


.66415 


0.88156 


1 


.6250 


6.9799 


4. 


2953 





.045139 


2 


.66467 


0.88618 


1 


.7500 


7.4359 


4. 


,2491 





.041551 


2 


.66585 


0.89672 


1 


.8750 


7.9127 


4. 


.2201 





.038262 


2 


.66709 


0.90918 


2 


.0000 


8.4036 


4. 


,2018 





.035261 


2 


.66835 


0.92183 


2 


.5000 


10.447 


4. 


.1786 





.026166 


2 


.67304 


0.97944 


3 


.0000 


12.547 


4. 


.1822 





.020291 


2 


.67677 


1.0401 



TABLE IX: Solution sequence for the EOS FPS. 



surface, and during the inspiral phase of binary neutron 
star, it is considered to be conserved. The rest mass of 
one component of a binary system is written Mq and 
defined by 

M := / pu a dS a = I fntcnlPyfitfx (CI) 

where dS a = V a ty/—gd 3 x, and y/—gd 3 x = aip 6 ^f^jd 3 x 
— aip 6 r 2 sin 0drd8d<fr, because 7 = / is assumed. 

In this paper, the mass Mi is used to specify an equal 
mass BNS sequence, and M = 2M\ is used to normalize 



d/R 2d/M Rq/M SIM Madm J/M 2 
1.4375 7.8579 5.4664 0.038232 2.66765 0.91610 
1.4062 7.7372 5.5020 0.038958 2.66749 0.91413 
1.4688 7.9817 5.4344 0.037499 2.66786 0.91847 
1.5000 8.1086 5.4058 0.036772 2.66810 0.92130 
1.5625 8.3727 5.3585 0.035291 2.66867 0.92711 
1.6875 8.9226 5.2875 0.032465 2.66983 0.94095 
1.7500 9.2099 5.2628 0.031101 2.67046 0.94804 
1.8750 9.7970 5.2250 0.028611 2.67166 0.96421 
2.0000 10.402 5.2008 0.026341 2.67286 0.98037 
2.5000 12.914 5.1657 0.019492 2.67718 1.0507 
3.0000 15.493 5.1643 0.015068 2.68040 1.1206 

TABLE X: Solution sequence for the EOS BGN1H1. 



d/Ro 2d/M Rq/M SIM Madm J/M 2 

1.4375 6.0667 4.2203 0.053901 2.66264 0.86655 

1.4688 6.1649 4.1974 0.052878 2.66279 0.86785 

1.5000 6.2628 4.1752 0.051867 2.66297 0.86944 

1.6250 6.6732 4.1066 0.047891 2.66396 0.87718 

1.6875 6.8881 4.0818 0.045973 2.66451 0.88202 

1.7500 7.1090 4.0623 0.044109 2.66512 0.88709 

1.8125 7.3353 4.0470 0.042335 2.66576 0.89298 

1.8750 7.5659 4.0352 0.040622 2.66635 0.89878 

2.0000 8.0365 4.0183 0.037449 2.66765 0.91094 

2.5000 9.9949 3.9980 0.027806 2.67237 0.96616 

3.0000 12.009 4.0031 0.021559 2.67620 1.0244 

TABLE XI: Solution sequence for the EOS ALF3. 

quantities. M\ is the gravitational mass of a single spher- 
ical star whose rest mass is equal to the rest mass Mq of 
one neutron star in the binary system of each model (see 
Table HH. 

The ADM mass Madm is rewritten using conformal 
spatial metric, 

Madm := ^ J (f ac f bd -f ab f cd )D b J cd dS a 

'■= 7^— / {f ac f bd -.f ab f cd )D b lc d dS a 

+ -L / (-2)f ah D b ^dS a 

= -7T" / D a i)dS a , (C2) 

where, in the second equality, the first term vanishes be- 
cause of our choice 7 = /; and ip — > 1 is used in the 
second term. We have calculated approximate values of 
Madm using this surface integral at the boundary of the 
computational domain. Also, we fit M^/2r to ip — 1 near 
the boundary, to ensure a constant Mw, ~ Madm ■ In the 
tables, however, the values of Madm are calculated from 
a formula in which the above surface integral is converted 
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to a volume integral using the Gauss-Stokes lemma. We 
apply this on the conformal spatial hypersurface, which 
results in a simpler formula; since, at spatial infinity ip — > 
1 ; ^ab ^ /0 6 and dSa = v a r^fd 2 x = V a r^d 2 x =: dS a , 
we have 



mass [40(, Madm = Mk- The equality is related to the 
relativistic virial relation for the equilibrium 41| , 



(C6) 



DM 



-!- / D a ipdS a , = --!- / AipdS, 



2n 



2n 



2n 



2mp 5 p H 



(C3) 



In the WL/NHS formulation the asymptotic fall-off of 
each field is sufficiently fast to enforce the equality. And 
in this case, the above two definitions for Mk agree as 
well. 

Finally, the total angular momentum is calculated from 
a volume form of surface integral at spatial infinity 



The Komar mass associated with a timelike Killing 
field t a is written 



~ / K%ct> b dS a = -!- / K\4> b dS a . (C7) 



1 

An 



-— I y a t p dS, 



a/3 



— I R%t P dS, 



An 



(2T a p-Tg a p )tf } dS a , 



J := 



To calculate J, we set the surface near the boundary of 
the computational domain of the central coordinates and 
use the Gauss-Stokes lemma to write 



:{pK + S)-2 ]a (3 a ]^^d? Xl (C4) J = ^ J^D a {K\4> b )dS 



where dS a — n ay /jd 3 x is used. To derive this, the global 
existence of a timelike Killing field is assumed. For the 
spacetime of WL/NHS formulations, no such timelike 
Killing field exists. Instead, an asymptotic Komar mass 
can be written 



1 
1 

8^ 



(87rj a r + K\D a ^ b ) dS. 



87Tj a a + A a b D a 4> b + —K(j) a D a tp 

'i[j 6 ^d 3 x. (C8) 



My, 



An 



Wt dS, 



a/3 



1 

-in 



D a adS a 



— I An <1}2 
An 



1 

A^. !i: 



aA ab A ab + -K 2 



+ Ana (/9 H + S) 



(C5) 



where (G a p — SnT a p)g al3 = is used. 

In [l0l | . we have derived asymptotic conditions for an 
equality of the ADM mass, and the asymptotic Komar 



The values of J listed in the tables in next section, are 
calculated from the latter formula. 



APPENDIX D: SELECTED SOLUTION 
SEQUENCES 

Selected waveless solutions of irrotational BNS for 
parametrized EOS presented in Table IIIII of Sec. IIVI are 
tabulated. All quantities are dimensionless in the geo- 
metric units G — c — 1 , except for the ADM mass which 
is in a unit of solar mass Madm [M q ]. 
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